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CHAPTER 1 


INTRODUCTION 
1.1 The Microburst 

The term downburst was used by T. Theodore Fujita to describe a form of low- 
altitude windshear in which a downdraft of air impacts the ground and diverges 
into a horizontal outflow. Outflows with a radius greater than 4 km are called mac- 
robursts while events with smaller outflows are called microbursts [1]. Low-altitude 
microbursts have been recognized as a hazard to aviation for aircrafts attempting to 
fly through this event. An aircraft encountering a microburst will first experience 
a performance increasing headwind as it enters the forward outflow, but encounters 
a strong downdraft as it approaches the center of the microburst. The performance 
decrease caused by the downdraft is further increased by the tailwind created as the 
aircraft flies out of the microburst. In some cases the performance decrease caused by 
such an event would render an escape impossible. A microburst event is represented 
in Figure 1.1. 

1.2 Windspeed Gradient 

Typical procedures for airborne windshear detection rely on the estimation of the 
windfield characteristics present within a region of space positioned in front of the 
aircraft. This space represents the region of possible flight paths for the aircraft. A 
potentially hazardous condition exists in the forward-looking direction when there is 
a large headwind-to-tailwind gradient in the windspeed. Such a gradient typically re- 
flects the windspeed characteristics observed by an aircraft as it penetrates the outflow 
of a microburst event. A sampling of average windspeed measurements throughout 
the outflow region will typically follow a gradient in windspeed feature known as an 
“S”-curve pattern. An example “S”-curve is shown in Figure 1.2. Note the rapid 




Figure 1.1 Symmetric Microburst 
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change from a negative windspeed (headwind), to a positive vvindspeed (tailwind). 
From these measured windspeed values, an attempt can be made to assign a value to 
the effect of the corresponding windfield conditions upon the performance level of the 
aircraft as it flies through the given region of space. One such measure of a hazardous 
situation has been defined in [2] in terms of a maximum windspeed divergence of over 
10 meters per second (m/s) occurring in a region of less than 4 km. 

It is important to note that the windfield estimates from forward-looking airborne 
sensors, such as pulse Doppler radar, LIDAR (Laser Interferometric Detection and 
Ranging) [3], and FLIR (forward looking infrared radiometer) [4] have only the ability 
to estimate radial (line-of- sight) windspeeds along the direction where the sensors are 
pointed. Current research continues on the ability of these systems to provide a 
reliable advanced warning of a windshear condition. 


1.3 Hazard Factor 

It is necessary to define a measure of the severity of a microburst in terms of its 
performance-decreasing effect on an aircraft flying through the event. The analysis 
involved in computing such an index relies upon the concept of airplane total energy 
which is the sum of the aircraft’s kinetic and potential energy [3]. Temporally and 
spatially varying windfields will contribute to the rate of change of airplane total 
energy. Severe windfield events, such as a microburst, can significantly lower total 
energy to the point where the aircraft can no longer remain in the air. 

Define the airplane total energy as 

E = ^mi/ 2 + hmg (1.1) 


where m is the mass of the aircraft, v is airspeed, h is altitude, and g is acceleration 
due to gravity. The potential altitude is defined as energy per unit weight and is 
given by 


hp 


E v 2 , 
— = - — V h 
mg 2g 


(1.2) 
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Computing the potential rate of climb assuming energy loss is negligible yields 


E vv ■ 

= f- h 

mg g 


( 1 - 3 ) 


which, when combined with aircraft equations of motion results in 


■ h _ E_ _ T-D 

mg mg 


W x W h . 

cos 7 H sin 7 

.9 9 


Wh 

v 


(1.4) 


where T is aircraft thrust, D is drag, 7 is the flight-path angle, and W x and Wh 
are the horizontal and vertical components of the wind velocity. Assuming nearly 
stable flight, i.e., T « D and 7 % 0, equation (1.4) reduces to the hazard factor or 
“F-factor” equation 

F = Wh. 

9 v 

Equation (1.4) now becomes 


(1.5) 


^ = A = L_£_f„. (1.6) 

mg mg 

Equation (1.6) defines quantitatively the performance decrease that a windshear con- 
dition may inflict upon an aircraft. Note that the “F-factor” in this equation is 
derived solely from windfield information, providing an aircraft-independent measure 
of the severity of the windfield. The thrust-to-weight ratio for any type of aircraft is 
then factored in to provide an aircraft specific hazard index. One common measure 
of hazard is defined for “F-factor” values above a typical threshold of 0.15. 


1.4 Hazard Factor Algorithm 

A method is presented in [5] for computing a hazard factor from windspeed esti- 
mates based on a weighted least-squares technique that reduces the effect of noise due 
to stationary and moving ground clutter. This hazard factor algorithm is currently 
implemented in the Airborne Windshear Doppler Radar Simulation (AWDRS) pro- 
gram [6] in addition to all hazard factor computations presented in this thesis. The 
technique estimates the slope of a linear function of windspeed vs. range for a given 
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range. Each windspeed measurement is weighted by the spectral width of the en- 
ergy return, providing an estimate of the confidence of the measurement. Currently, 
windspeed estimates from five adjacent range cells are used for each computation of 
hazard. An advantage to the weighted least-squares algorithm is the addition of a 
measure giving the error of the estimate. This error is the averaged sum of the least- 
square residuals and can be used to identify wide cell-to-cell variations in windspeed 
estimates. Consider the measured windspeed, Vj over m range cells (m ~ 5), 

Vj = bi + b 2 r 3 + &j (1.7) 

where b x and bi are constants, Vj is the estimated mean Doppler velocity at range 
rj, and ej is the least-squares approximation error. The estimate 62 represents the 
velocity- to- range slope and is desired for the estimate of hazard. Begin the least- 
squares solution for 62 by representing (1.7) in matrix form 

v — Ab + e (1.8) 


where v is an m x 1 column vector containing values of vj, b is a 2 x 1 column vector 
containing b\ and hi respectively, e is the column vector of values, and A is an 
m x 2 matrix of the form 

1 r x 
1 r 2 


A = 


(1.9) 


By adding a weighting matrix to the least-squares solution of 1.8, the estimates of Vj 
with the greatest confidence level will influence the least-squares solution more, and 
the estimates which have a low confidence level will not influence the solution of (1.8) 
as much. The weighting matrix $ is composed of spectral width measurements cr^, 
producing 



(1.10) 
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A discussion of spectral width estimators follows in Chapter 2. The weighting matrix 
is included into the least-squares pseudoinverse solution of (1.8) for matrix b [7] 


6= (a t $ _1 a) 1 A t $- 1 u . 


(1.11) 


From (1.11) the hazard level is contained in the b 2 component of matrix 6. Smoothing 
of this estimate over adjacent range cells through averaging is possible, although 
results presented in this thesis are without averaging. The error estimate or “goodness 
of fit” is equal to the weighted sum of the error vector e 


r = 



m 


( 1 . 12 ) 


and can be used to exclude an estimate as invalid if r is greater than a given threshold, 
indicating that the least-squares fit has substantial error. 


1.5 Doppler Weather Radar 

The advantage of pulse Doppler radar over conventional radar is Doppler’s ability 
to identify not only range to an object in the radar’s scanning area but also the relative 
velocities of these objects. For this reason, Doppler weather radar may prove to be an 
effective device in the detection and avoidance of hazardous low-altitude windshear by 
identifying approaching windfield characteristics. The ranging capability of the radar 
is achieved by measuring the time in which a transmitted pulse of energy is reflected 
off objects and returned to the receiver. The operating characteristics of Doppler 
radar from a time domain perspective have been thoroughly presented [8, 9, 10]. 
Here a frequency spectrum analysis of Doppler Radar is provided. 

Scatterers which return transmitted radar energy will induce a Doppler shift 
in the frequency of the transmitted pulse energy proportional to the radial velocity 
between the object and the transmitting/receiving antenna. Objects approaching the 
antenna will cause a Doppler shift of the transmitted energy to a higher frequency 
while a diverging object relative to the antenna will create a Doppler shift to a lower 
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frequency. Consider a transmitted pulse. s(t), having amplitude ,4 which is composed 
of a single frequency, / 0 , 

s(t) = Acos(2irf 0 t) (1.13) 


which has the following frequency spectrum: 

S(f) = F[s{t ) ] = \-Ae^ fot+ ^ + -Ae- j{27rfot+ * ] . (1.14) 


The amount of transmitted energy returned to the receiving antenna is equal to the 
sum of the return energy from all scatterers. The frequency of the return from scat- 
terer i, will be Doppler shifted by an amount proportional to the relative velocity 
between the receiver and the scatterer and the total Doppler spectrum will be a 
distribution of returns from all scatterers. The relationship between the radial veloc- 
ities of the scatterers and the corresponding shift in Doppler frequency can be shown 
through the Doppler equation for light as predicted by the theory of relativity [11] 



where / is the transmitted frequency and f is the frequency observed by an object 
where the relative velocity between the source and the target is u and c is the speed 
of light. Negative u indicates that the objects are converging and positive u indicates 
that the objects are diverging. In the case of Doppler radar, the Doppler shift must 
be computed twice: first when the transmitted pulse reaches a scatterer, and again 
when the returned energy reaches the radar antenna. At the scatterer, the Doppler 
shift of the transmitted pulse can be described by 


1 _ t'i-Hq 

c 

1 -(^) 2 

where f- is the frequency observed at scatterer i, and u, and v a are the radial velocities 
of scatterer i and the aircraft, respectively. The aircraft velocity in included in order 
to remove the effect of a moving radar platform. The next step is to compute the 



(1.16) 
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Equation (1.19) describes the relationship between the radial velocities of an airborne 
Doppler radar system and scatterers in its range cell volume and the Doppler shift 
in frequency of the returned signal. Assuming a transmitted pulse having the energy 
spectrum given in (1.14), the spectrum of the returned energy is 

Sret(f) = £ ( 1 . 20 ) 

where A , is the amount of energy returned from scatterer i and fa is the phase offset 
due to scatterer i. Because the returned energy contains Doppler shifts about ±f 0 , 
it is necessary to demodulate the transmitted signal in order to obtain Doppler fre- 
quency shifts that are centered about zero frequency, as in stationary objects, rather 
than the transmitted frequency f 0 . This demodulation is accomplished by multiply- 
ing the returned signal, S re t by the single-sided complex exponential at frequency 
f a — f 0 (l 4- 2u a /c), resulting in a negative shift in all frequencies of S ret by f s . Note 
that the modulation term f 3 includes the aircraft velocity, which must be known. The 
shifted signal becomes 


SlQ(f) = SrM)*SU-f.) 




( 1 . 21 ) 

( 1 . 22 ) 
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Applying (1.19) and the formula for the wavelength of an electromagnetic signal. 
A 0 = c/f ot to (1.22) yields 

Siq(f) = 2 + ^A t e- J (M 2 A- %]'+*•) } . (1.23) 

The second term in (1.23) can be removed through a low-pass filter resulting in 
a complex baseband signal representing the returned Doppler shift centered about 
zero frequency. This signal is commonly represented by its resulting in-phase and 
quadrature-phase time series returns of the form 


1 - [ Ai cos + &)] ( L24 ) 

Q = [ A » siti ( -p - 1 * + </>»•)] • (1-25) 

Equations (1.24) and (1.25) define the complex time series referred to as IQ data 
in terms of the sum of Doppler shifts caused by all scatterers. I refers to the real 
component of the complex series while Q refers to the imaginary component. For 
any single scatterer, the complex frequency of the IQ data will be related to the 
line-of-sight velocity of the scatterer relative to ground, that is, 

fdo PP = -y (1.26) 

where v is the ground velocity of the scatterer projected along the line-of-sight from 
the radar antenna to the scatterer. For a given sampling rate equal to the pulse- 
repetition frequency (PRF), the greatest unaliased frequency returned is 


f — —— 

J max 


(1.27) 


and the maximum unambiguous velocity of a scatterer is 

A 0 


Umax — 


iT, 


( 1 . 28 ) 


The radar parameters used for all analysis include a PRF of 3755 pulses-per-second, 
a pulse width of 0.96/ns, and a transmitted frequency of 9.33625 GHz (A„ = 3.21103 
cm). From (1.28) the maximum unambiguous relative velocity for this situation will 
be equal to 30.14 m/s. 
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1.6 Problem Statement 


The ability to detect a hazardous condition with the ‘‘F-factor" index relies 
heavily upon the accuracy of the windspeed estimates obtained. Strong clutter returns 
can heavily bias pulse-pair and spectral mean estimates of windspeed. A common 
approach to increasing the signal-to-clutter ratio for better windspeed estimation is to 
implement a notch filter to reduce the main clutter return which is usually present near 
zero relative velocity. This filtering can, however, attenuate the weather return which 
will further reduce the accuracy of the windspeed estimate. Proposed is a method 
whereby an unbiased estimate of windspeed can be achieved in a computationally 
attractive procedure and without clutter rejection filtering. 
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CHAPTER 2 


WINDSPEED ESTIMATION 

Although a Doppler weather radar system cannot directly measure windfield con- 
ditions, it is possible to infer information about such conditions when the radar return 
is a backscattering of energy from targets which are suspended particles and whose 
motion is controlled by the prevailing windfield conditions. As shown in section 1.5, 
the returned radar energy for a given volume will consist of a collection of returns 
from all targets that scatter or reflect the incident radar energy. In the case of a 
Doppler weather radar system, the targets are assumed to be particles suspended 
in the resolution volume, such as rain droplets, dust, or insects. From the Doppler 
spectrum, which contains a distribution of returns from all targets, windspeed estima- 
tion techniques attempt to characterize the windfield motion in the current resolution 
volume based on the motion of the targets in the resolution volume. 

2.1 Time Domain Pulse-Pair Mean and Width Estimator 
The time domain pulse-pair processor is a mean frequency estimator that pro- 
vides an estimate of Doppler spectral mean and width in terms of the first autocor- 
relation lag. The pulse-pair Doppler mean estimate is [12] 

( 2 . 1 ) 

with the autocorrelation of the time series at lag T s given by 

1 N ~ 2 

R(T t ) = -rr Y, x*(n)x(n 4- 1) . (2.2) 

jV n= 0 

Let x(n) be a series of N complex IQ data samples recorded at the pulse repetition 
interval, T s , where T s = 1 / P RF . Equation (2.2) is then equal to the complex auto- 
correlation of x(n) at a lag of one pulse, thus, define Ri = R(T S ). The corresponding 



pulse-pair Doppler width estimate for x(n) is 


(Tv = 


In 


' S 

J«il 


1/2 


sgn 


2kT s V2 

In (2.3), S is the total power in the data sequence, 

JV n=0 


In 


}Rx 


(2.3) 


(2.4) 


The sgn term is included in order to tag the width estimate with a negative value on 
the occasion when the estimate R(T S ) is greater than the estimate S. For the case 
where < 1, an alternate width estimator [12] can be used 

d/2 / 


<t„ = 


2 nT s y/2 


1 - 


|i?i 


sgn 1 - 


l^i 


(2.5) 


This estimate of spectral width is only valid over the given range of narrow widths 
due to a bias present at larger spectral widths. 


2.2 Frequency Domain Pulse-Pair Mean and Width Estimator 
The frequency domain pulse-pair mean estimator performs a similar analysis with 
similar results as its time domain counterpart. The difference is that the frequency 
domain estimator computes the autocorrelation of the sequence at lag T s from the 
power spectrum rather than from the time series. The equation for Doppler mean 
estimate remains the same as in (2.1), that is, 

v = ^r*.rgR(T,). (2.6) 


However, the first autocorrelation lag can be derived directly from the power spectrum 
with 




N 


£ W{k)e j2 7r 


k=0 


(2.7) 


where the power spectrum is defined as 


»'(<■)= |.V(A-)| 2 = 


A T — 1 

E 

71=0 


x(n)e J a ' 


( 2 . 8 ) 
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In order for the autocorrelation of (2.7) to yield the same results as (2.1 ) at least one 
zero must be appended to the time series x(n ) before the discrete fourier transform 
(DFT) is implemented to give X(k). This is due to the nature of the DFT, which 
assumes a periodic time series. Without a padding of at least one zero value, the 
inverse DFT implementation of (2.7) would correlate the last sample x(N — 1) with 
the first sample a-(O). The frequency domain spectral width estimator resembles the 
corresponding time domain width estimator given in (2.3) and (2.5) except for the 
modified autocorrelation estimate R i and the total power estimate S. The total power 
is equal to the sum of the power spectrum 

(2.9) 

iV k=0 

The computation of pulse-pair spectral mean and width in either the time domain 
or the spectral domain yield similar results. Therefore, the implementation of one 
technique over the other would normally be based on computational costs. 


2.3 Spectral Domain Mean-Variance Mean and Width Estimator 
An alternative algorithm estimates the Doppler mean and width by computing 
the mean and variance of the given power spectrum directly in the spectral domain. 
For this algorithm, the first and second moments about an estimate mean bin number, 
k m , are computed [12]. k m is determined from the frequency bin in W(k) containing 
maximum energy. The Doppler mean estimate is then computed with 


v = 


-A 

2A r 


km + 4=r E (* ~ M W[mod N (k)} 


ST. 


Sk=km~i 


( 2 . 10 ) 


In (2.10), S is the total power of the data sequence and mod s(k) denotes k modulo 
N. The Doppler width estimate a v is equal to the variance of the power spectrum 
about k m , 


cr,. = 



1V2 


W[mod.v(fc)] 


( 2 . 11 ) 
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2.4 The Effect of Clutter on Windspeed Estimates 

A major concern with airborne pulse Doppler radar systems is the presence 
of strong energy returns from targets not controlled by wind. These targets can be 
stationary or moving objects on the ground, or objects in the range cell volume whose 
motion is not controlled by the wind such as birds or other aircraft. In situations 
involving heavy rainfall, it has been observed that the main rain shaft may not be 
under the control of divergent winds and may appear as a mode about Doppler zero, 
similar to a ground clutter signature. Unwanted returns originating from the ground 
are called ground clutter and can significantly inhibit a system's ability to correctly 
identify a hazardous condition. Due to physical requirements which necessitate the 
use of antenna elevation angles directed at or near the ground, airborne weather 
radar systems are more susceptible to the effects of ground clutter contamination 
than ground based systems with much higher antenna elevations. 

Although microburst conditions can occur at any altitude, it is the low-altitude 
windshear events which inherently pose the greatest hazard to an aircraft penetrating 
such an event due to the horizontal shear created by the microburst as it impacts the 
ground. This hazard is also complicated by the low altitude of the scenario which 
provides the pilot with a relatively short response time in which to avoid catastrophe. 
In order for an airborne radar system to observe such low-altitude events, antenna 
elevation angles must be kept small so that any hazardous windshear conditions can 
be identified. Antenna elevations currently being used are typically within three 
degrees above or below horizontal. At these elevations, radar return from antenna 
sidelobes can introduce much unwanted clutter into the received signal, especially at 
closer ranges. 

Moving objects on the ground such as cars on a nearby interstate highway can 
create "discrete clutter which will appear as a narrow spectral peak within the re- 
turned IQ data. This form of clutter can impair the ability to correctly estimate 
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windspeed values because the clutter can be located anywhere in the Doppler spec- 
trum and thus, is difficult to characterize. The implications of discrete clutter on 
modal analysis are discussed in Chapter 4. 

Because clutter from stationary reflectors can be characterized as a spectral mode 
near Doppler zero, it can be more readily removed through a clutter rejection filter. 
The presence of a strong clutter return can severely impair the ability of mean spectral 
estimators such as pulse-pair and mean-variance to correctly identify the windspeed 
conditions present in a given range cell. Clutter can bias a windspeed estimate toward 
Doppler zero or in the event of a strong clutter return and a weak weather return, 
completely cover the weather return signature. However, by identifying clutter char- 
acteristics, measures can be taken to reduce the unwanted effects of clutter upon 
windspeed estimators. 


2.5 Clutter Rejection Filters 

The general approach to reducing the effect of clutter on windspeed estimates 
is to remove the clutter return through filtering. Several types of clutter rejection 
filters exist which can improve the quality of mean spectral estimators by attempting 
to remove returns from clutter targets while retaining returns from weather targets. 
Because weather and clutter returns may overlap both spatially and temporally, a 
complete removal of clutter return energy without some attenuation to any weather 
return present is not possible. In general, clutter energy is characterized as a nar- 
row band near Doppler zero. Clutter rejection filters are typically notch filters which 
remove spectral energy near Doppler zero. Some typical clutter rejection filtering 
methods include (1) finite impulse response (FIR) notch filter [13. 14. 15]. (2) infinite 
impulse response (HR) notch filter [13. 14. 15]. and (3) spectral domain line editing 
(ideal notch filtering) [13]. The disadvantages to clutter rejection filtering are the ad- 
ditional computational resources needed and the possible attenuation of the weather 
signal below a detectable level. 
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CHAPTER 3 


SPECTRAL ESTIMATION USING MODAL ANALYSIS 

Typically, mean estimators such as the pulse-pair method function by attempting 
to reduce the energy return through filtering so that the remaining return appears to 
have originated from a single source. 

Proposed is a method for windspeed estimation based on a modal analysis ap- 
proach which would not require filtering of an unwanted clutter return. Unlike the 
pulse-pair estimator, a modal analysis approach attempts to model the returns from 
all sources. In this way, it attempts to identify all dominant sources of energy return 
present without having to rely on a single mean estimator for the windspeed such as 
the pulse-pair estimator, which can only yield a single spectral location parameter 
associated with the return data. When the spectral content of the return data in- 
cludes more than one mode, the single pulse-pair estimate will attempt to reflect both 
modes with a single value producing an estimate that is potentially uncharacteristic 
of either mode. To compensate for this, clutter rejection filtering generally is used in 
conjunction with the pulse pair-estimator. Filtering attempts to remove the return 
energy from undesired sources while retaining the return energy from the desired 
weather source. In contrast, a modal approach to windspeed estimation models the 
return energy from all sources and then attempts to classify the source of each return. 

For the problem at hand, the primary sources of return energy to the Doppler 
radar system are weather and clutter. Based on the assumption that only these 
two primary sources of return energy exist, a method of spectral estimation which 
involves the modeling of two dominant modes is sufficient. Once the identification of 
the spectral modes is achieved, a statistical two-class pattern recognition routine can 
be employed for the classification of modes. 



3.1 The Prony Approach 


The Prony method was created by Gaspard Riche, Baron de Prony as a means 
of representing the expansion of various gases in terms of the sums of complex ex- 
ponentials. Assume that there exists a sequence, x n , consisting of N complex data 
samples. The Prony method generates an alternate series, x n , which approximates 
the original series through a summation of M complex exponential terms [16, 17]. For 
an approximation where N < 2 M, the exact series can be represented and x n will 
equal x n . For situations where N > 2 M, an approximation of x n containing some 
error will be generated. Define the approximation 

M 

in = CmHn, for n = 0, • • ■ , A r - 1 (3.1) 

771 = 1 

where C m and are complex values of the form 
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so that the values of p m in (3.1) correspond to the roots of equation (3.6). Next, 
multiply (3.1) by a p and substitute (n — p ) for n so that 

M 

x n - r a„^'£C m a l ^ m -<’ (3.7) 

m=l 

for 0 < n — p < N — l. Next, replace fi^ p with to get 

M 

%n-p a p = ^2 Cm Up V • (3*8) 

m= 1 

After the summation of both sides of (3.8) over 0 < p < M note that the final 
summation of 


M 

y! %n—p &p 
p=0 


M M 


m~l 


p=o 


(3.9) 


is of the form of (3.6) evaluated at one of its roots and thus, equals zero. The left 
side of (3.9) will form the recursive difference equation 

M 

X n — ^ '] dm Xn—m (3.10) 

771=1 

defined over m < n < N — 1. Note that (3.10) is in terms of the approximation series 
x which is still undetermined. By including the approximation error e from (3.5), 
the difference equation of (3.10) can be expressed in terms of the original series x n . 
Substituting for the approximation series in (3.10) results in 

M 

x n — ^ d m x n — m -f* e n (3.11) 

m= 1 

M M 

= &m x n—m “h } ^ &m €-n—m (3.12) 

m=l 77i=0 


defined for M < n < N — 1. The Prony method requires the solution of (3.12) for 
the a’s while minimizing the squared error e n . An iterative solution to this difficult 
nonlinear system has been presented by Huggins and McDonough [18]. 
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3.2 Extended Prony Analysis 


The extended Prony technique is a method whereby the nonlinear characteris- 
tics of (3.1) can be linearized at the cost of an acceptable, yet sub-optimal solution 
[16]. The error measure for the Prony technique has been shown to be the difference 
between the sampled data x n and a linear prediction series based on M previous 
samples of the approximation series. To simplify (3.11) into a linear equation where 
the values a m can be determined efficiently, define a modified error measure equal to 
the difference between x n and a linear prediction series based on M previous samples 

%n—M+ 1 1 • • • 1 %n—l } 5 

M 

X n = ) d m %n—m T (3.13) 

m= 1 


where the modified error is 


M 

£ n = ^ ^ &n—m • (3.14) 

m—0 

By minimizing this revised error measure, (3.13) degenerates into an efficient spec- 
tral estimation technique known as Autoregressive (AR) parameter estimation. Ap- 
pendix A describes AR modeling and presents an efficient method for AR parameter 
estimation. 


3.3 Second Order Extended Prony Analysis 
Presented in this section is an implementation of a second order extended Prony 
technique applied to windspeed estimation of pulse Doppler weather radar data. The 
extended Prony approach begins with any AR parameter estimation technique which 
will generate the second order coefficients a i and a 2 . The Levinson-Durbin algorithm 
described in Appendix A will be used to generate the AR coefficients due to its low 
computational requirements which would facilitate implementation on a real time 
windshear detection system. From the AR coefficients a i and o 2 the values of fi\ and 
// 2 from (3.1) are determined through rooting the characteristic polynomial 

+ d\/.i -f ci 2 = 0 . (3.15) 
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Note that pi and fi 2 are the pole locations in the z-plane of the AR process. For the 
second order case, these pole locations can be solved by way of the quadratic equation 


P+ 

P- 


—di + \Ja\ - 4a 2 
2 

-aj - \Ja\- 4 a 2 

_ 


(3.16) 

(3.17) 


The desired frequency estimates as well as damping factors are derived from 
the values of fi+ and /i_. The parameters C\ and C 2 contain information about 
amplitude and phase of the Prony model, but are not required for an extended Prony 
analysis. Recalling that the p’s are related to the frequency estimates by p m = 
exp [{a m + j2irf m )At] and that At is equal to the inter-pulse period T s , the frequency 
values fi and f 2 are computed with 

fi = arg (/ z ') • ( 3 - 18 ) 


The desired windspeed estimated can now be inferred from the conversion from fre- 
quency to velocity 

v i = ar g(/*«) • (3-19) 

Note the similarity between (3.19) and the pulse-pair equation given in (2.1). In fact, 
the pulse- pair estimate of Doppler velocity is analogous to a windspeed estimate from 
a first order AR model where the argument term would be the coefficient a]. 


For reasons to be made clear in Chapter 4. a different assignment of the two 
poles will be defined according to the magnitude of the respective Doppler velocities. 


That is for pole 1 


and for pole 2 


f p+, Kl > K’-| 

( /<_, else 


P2 


/'+• l r -l > K’+l 

p_. else 


(3.20) 


(3.21) 
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The damping factors <*i and q 2 yield information about the spread of a spectral 
mode. Small damping factors indicate wide spectral modes while larger damping 
factors indicate narrower modes. The damping factors are equal to 


I/hi 

PRF ‘ 


(3.22) 


It has been stated that the extended Prony method can yield a greatly simplified, 
yet sub-optimal solution compared to the Prony technique. The concession made 
by the extended Prony method is that the damping factor, a TO , and amplitude, A m , 
have been reduced to a single parameter. This parameter is directly related to both 
the damping factor and the amplitude of the spectral mode. A pole value with a 
magnitude near unity will have a relatively large magnitude and a relatively small 
damping factor, while a pole with a small magnitude will produce a smaller amplitude 
and a larger damping factor. When one factor is changed the other will be changed 
as well. This concession is not a concern in the present windshear detection problem 
as only amplitude, and not the damping factor, is used in the current analysis. 

A comparison between a second order AR model spectral estimate and an FFT 
spectral estimate for a sample data set is provided in Figure 3.1. This 96 point 
IQ data set contains returns from two distinct spectral modes: one mode is centered 
about zero Doppler velocity while the other mode is at some negative velocity indi- 
cating motion in a direction toward the aircraft. For the FFT operation, the sample 
data has been zero-padded to a length of 128. Pole 1 and Pole 2 reflect the values of 
Hi and H 2 converted into a polar form where the pole’s magnitude and argument are 
represented. Note how the magnitude of each pole is reflected in the relative height 
of the spectral peaks and that each pole’s argument converted to Doppler velocity 
through (3.19) corresponds to the Doppler velocity of the spectral modes and the 
local maxima in the AR spectral density. Even though a second order model such 
as the the example in Figure 3.1) cannot be expected to provide an accurate and 
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Figure 3.1 Second Order AR Spectral Estimate and Pole Locations 



detailed spectral density, it is possible, however, for such a model to provide accurate 
estimates of the frequencies associated with dominant spectral modes. 

It is now possible to combine the pole values from a single range cell, such as 
in Figure 3.1 with other pole values from all range cells in any given radar snapshot. 
Figure 3.2 is a bubble plot which represents a series of 91 range cells and the pole 
values which correspond to that cell. From the radar parameters given in section 1.5, 
each range cell represents approximately 145 meters. On the bubble plot each range 
cell has two circles associated with it which represents the argument and magnitude 
of the second order pole locations. Pole magnitudes, that is their distance from the 
xr-domain origin, is shown by the size of the circles, while the argument of a pole is 
translated into its corresponding Doppler velocity estimate as computed by (3.19). 
When viewing the two velocity estimates associated with a second order extended 
Prony technique for the data in Figure 3.2, two distince modes appear indicating that 
the energy in the return is from two or more separate return sources. In Figure 3.2 
both a weather mode and a clutter mode are apparent. The “S”-curve feature can 
be classified as weather return from a microburst event [19] while the strong mode 
present throughout all range cells near Doppler zero can be classified as a clutter 
mode. In this example, range cells 65-96 contain both a strong central clutter mode 
as well as a weaker model located near ±30 m/s. The pole values for these range cells 
can be attributed to neither a weather nor a clutter return, but is a result of the AR 
modeling process which will be discussed further in Chapter 4. 

The windspeed estimates of the weather mode can be considered to be an unbi- 
ased estimate because both spectral modes are modeled independently of each other. 
The clutter return should be reflected in the clutter mode, leaving the weather mode 
to reflect mostly weather return. For this reason, clutter rejection filtering is not 
required for windspeed estimation with this modal approach. 

A complication to this second order approach arises due to the generation of two 
velocity estimates for each range cell volume, but the need of only a single windspeed 
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estimate. Chapter 4 addresses this topic and presents various methods which attempt 
to resolve the problem of characterizing the windfield through a second order extended 
Prony approach. 
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CHAPTER 4 


WINDSPEED ESTIMATION FROM A TWO POLE MODEL 

As introduced in Chapter 3, a second order Prony model may avoid the need 
for typical clutter rejection filtering. However, further processing is required in order 
to formulate a valid windspeed estimate. The identification of two spectral modes 
by the extended Prony approach presented in Chapter 3 introduces the need for 
some form of classification which can distinguish those modes that identify clutter 
sources from those modes that identify weather sources. Presented in this chapter 
are two classification schemes which attempt to derive a windspeed estimate given 
the parameters of a second order Prony model. The first classification scheme forms 
a windspeed estimate based on the Doppler velocities associated with the poles of a 
second order all-pole model. The second scheme implements an approximation which 
allows for a windspeed estimate to be determined directly from the coefficients of a 
second order autoregressive (AR) model which reduces the computational load of the 
windspeed estimator. 

4.1 Characterizing the Poles 

Before attempting to classify spectral modes based on second order extended 
Prony method pole values, discriminating characteristics are needed which will exploit 
differences between clutter returns and weather returns. Such a statement assumes 
that dominant weather and clutter modes present in the Doppler spectra will be 
represented by a second order AR model. The complex pole values associated with 
such a model yield information about the Doppler velocity and relative intensity of the 
spectral modes through the pole s argument and magnitude, respectively. To examine 
how a second order model might represent weather and clutter characteristics, first 
analyze the Doppler spectrum over a series of range cells for situations with and 



without a weather return present. A 20th order AR model is used in Figures 4.1 and 
4.2 to provide a high order approximation of spectral density for each range cell in both 
a simulated dry microburst event (Figure 4.1 ) and a clutter only situation (Figure 4.2). 
Both situations simulate an aircraft on final approach to Denver Stapleton airport 
runway 26R with the microburst model D51. For each run the antenna azimuth was 
set at -10 degrees, antenna elevation was 1 degree above a glideslope of 3 degrees, and 
the aircraft was positioned 7.4 km from touchdown. All simulation parameters have 
been included in Appendix C. Data for the microburst model have been derived from 
data collected during an actual microburst event which occurred on July 11, 1988 
at Stapleton airport [20]. In the dry microburst situation, the energy return from 
the leading edge of the outflow can be seen as a negative Doppler shift(headwind) 
with a peak Doppler velocity of about —14 m/s in the vicinity of range cell 22, 
while the backside of the outflow has a weaker return but can be distinguished as a 
positive Doppler shift in the following range cells. In comparison, the “clutter only” 
situation, shown in Figure 4.2 does not contain the “S”-curve pattern commonly 
described as characterizing a microburst event. Instead, the return with little or 
no weather appears as a single spectral mode near Doppler zero. To examine the 
pole values for both the weather plus clutter situation and the clutter only situation, 
bubble plots similar to Figure 3.2 in Chapter 3 are included in Figures 4.3 and 4.4. 
It can be observed that for the dry microburst situation (Figure 4.3), the front edge 
of the microburst is reflected by negatively shifted poles in range cells 20-28 while 
the weaker outflow shows up as positively shifted poles in range cells 30-46. For the 
clutter only situation (Figure 4.4), there appears to be no pole pattern which would 
indicate any return other than clutter. Range cells 6-22 each contain one strong pole 
near Doppler zero and one weaker pole near ±30 m/s. These weaker poles (defined 
as pole 1 because they have the greater Doppler velocity) may at first appear to 
represent a weak weather mode since clutter is primarily near Doppler zero, but upon 
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inspection their presence can be attributed to a pole spreading effect resulting from 
the AR modeling process. 

4.2 Pole Spreading Effect from AR Modeling 

When the spectral content of a range cell is comprised primarily of a single strong 
clutter mode, the pole values computed with a second order extended Prony method 
will tend to double up to reflect this spectral mode. When the spectral content of 
the return is of a more even distribution, that is, no single spectral mode dominates 
the return, the poles will tend to spread out across the Doppler spectrum so that one 
pole is at the central clutter mode while the other pole is situated approximately half 
the distance around the unambiguous Doppler velocity spectrum. 

This spreading effect is detrimental to the classification of weather poles and 
thus, windspeed estimation. Difficulties arise due to similarities between weather 
poles and poles resulting from a modeling effect when no weather return is present. 
Much of the effort to classify weather poles from non-weather poles is directed toward 
distinguishing and removing poles caused by spreading from valid weather poles. 

The spreading out effect of the poles along the Doppler spectrum is not only 
limited to a second order model, but is also seen in higher order models as well. 
For example. Figure 4.5 is a bubble plot of the same clutter only simulation but 
with three poles per range cell. Observe that in range cells 6-22 the three pole 
values are distributed evenly throughout the Doppler spectrum, that is all poles are 
positioned approximately 20 m/s from the other two pole location in the Doppler 
velocity spectrum. These are the same range cells where pole spreading was observed 
in the two pole case. In the latter range cells, pole spreading occurs, but in this case, 
the central clutter mode is more dominant and attracts two pole values while only 
the third pole spreads away from the other two. 
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4.3 Two Class Pattern Recognition Technique 

In statistical pattern recognition terms, the second order extended Prony method 
sets up the basis for a two class pattern recognition problem [21] where the two classes 
of interest are weather and no weather. Since a windspeed estimate is actually the 
desired result, certain conditions will be incorporated into the classification system: 

1. Pole 1 is defined as that pole with the greatest Doppler shift. 

2. Pole 2 is defined as the remaining pole in each range cell. 

3. Pole 1 may or may not be classified as weather. 

4. Pole 2 will always be classified as non-weather. 

5. Both poles may be classified as non- weather poles. 

6. The windspeed estimate for a range cell with a pole classified as weather will 
be the Doppler velocity associated with the weather pole. 

7. The windspeed estimate for a range cell without a pole classified as weather 
will be zero. 

Condition 7 deserves further discussion. It states that if neither pole is given the 
classification of weather then the windspeed estimate will be assigned a zero value. 
For a return consisting of stationary clutter with pole spreading, or for a return with 
discrete clutter, a spectral mode away from zero will be present that is not a result 
of a weather return. A classification scene will attempt to distinguish between these 
non-weather pole sources and true weather sources in an effort to estimate windspeed. 
When pole 1 is away from Doppler zero and the classifier determines that it is not 
weather, then the windspeed estimate will be zero on the basis that no pole could be 
identified as a weather source. 

At this point, the definition from Chapter 3 that pole 1 is that pole with the 
greatest Doppler velocity shift begins to become apparent. The reasoning behind 
such a definition is due to the characteristics of weather and clutter modes. Because 
clutter modes are typically near Doppler zero and weather modes can be shifted some 
amount away from Doppler zero, the assumption can be made that if one of the 
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spectral modes is indeed a weather mode, then this mode will have a greater Doppler 
shift than the clutter mode present. Therefore, by classifying pole 1 and pole 2 in 
this manner, pole 2 will always have the classification of a clutter pole while pole 1 
may or may not be classified as a weather pole. 

Classification begins with the selection of features, or measurements used to 
distinguish between classes [21]. In this discussion, features will be limited to mea- 
surements of pole values independent of range cell number. Further research may 
incorporate pole values across a series of adjacent range cells for additional informa- 
tion which may enhance the ability of the classifier to distinguish between weather 
poles and non-weather poles. Taking this approach one step further, classification 
could draw upon previous classifications, which adds a time variable to the classifi- 
cation scheme. By increasing the scope of the classification process beyond a range 
cell-by-range cell decision, the negative effects of pole spreading and discrete clutter 
may be better eliminated. This extension to the pole classifier will be left as the topic 
of further investigation. 

The classification scheme presented in this discussion consists of a two-dimension- 
al feature space. The features chosen should be such that differences between the two 
decision classes are apparent, thus allowing for a decision to be made. One possible 
choice of features is the real component of each pole and the imaginary component 
of each pole. For a two-dimensional case such as the current estimation problem, 
plotting the two features against one another will usually give an indication as to the 
ability to form a reliable decision. Therefore, for the current problem, the feature 
space is the complex plane on which the values of pole 1 and pole 2 are to be plotted. 
Figures 4.6 and 4.7 are the values of pole 1 for both the simulated microburst event 
and the simulated clutter only event, while Figures 4.8 and 4.9 are the values of pole 2 
for the same situations. 

By comparing the values of pole 1 in Figures 4.6 and 4.7, any pattern differences 
between the two pole distributions may be attributed to the presence of a weather 
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Figure 4.7 Complex Values of Pole 1 for Simulated Clutter Only Return 
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Figure 4.8 
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Figure 4.9 Complex Values of Pole 2 for Simulated Clutter Only Return 
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return in Figure 4.6. Although a difference does exist in that the weather and clutter 
simulation produced poles with larger positive imaginary components than the clutter 
only simulation, an improved choice of features may exist. Since windspeed estima- 
tion, or the Doppler velocity shift caused by a weather mode is the final objective, 
choose a polar representation for the set of features that includes (1 ) Doppler velocity 
and (2) complex pole magnitude. By plotting pole values in this revised feature space 
the weather pole and the clutter pole attributes become more apparent allowing for 
the feature space to be partitioned into decision regions where classifications can be 
made. Figures 4.10 and 4.11 are distributions of pole 1 in the polar feature space and 
Figures 4.12 and 4.13 are distributions of pole 2 in the same feature space. 

Once a feature space is determined, the next step in classifying poles as weather or 
non-weather is to partition the feature space into decision regions where classifications 
can be made according to each pole’s location in the feature space. In comparing 
Figure 4.10 with Figure 4.11, a group of poles is present in the microburst simulation 
that is not present in the clutter only simulation. The poles are between —10 m/s 
and —15 m/s Doppler velocity with magnitudes greater than 0.6. Upon comparing 
the scatter plots in feature space with the bubble plots, it can be determined that 
the group of poles mentioned are the same poles which reflect the leading edge of the 
microburst event recorded in range cells 20-28. As seen in Figures 4.12 and 4.13, the 
Doppler velocities associated with pole 2 are in the region near Doppler zero. For 
this reason pole 2 is always given the classification of clutter. 

The decision regions which separate those poles classified as weather from those 
poles classified as non- weather will be determined through an analysis of pole distri- 
butions from a large collection of return data. The data used includes 10 simulated 
returns of the July 11th microburst event in Denver as well as 10 returns of the same 
situations without weather present. Although the decision regions determined from 
this data will be particular to the Denver area and the type of microburst included 
in the simulation, the approach to be discussed can be extended to any location or 
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Figure 4.11 Polar Values of Pole 1 for Simulated Clutter Only Return 
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weather condition. A more complete analysis may include data which is distributed 
among a variety of locations and weather conditions and may reveal the optimum 
decision regions given all situations. 

If weather poles and non-weather poles could always be separated due to their 
appropriate values then the classification scheme would be trivial. Unfortunately, 
weather poles and non-weather poles may occupy the same region in feature space. 
The decision boundary should be designed so as to minimize classification errors. 
Presented in Figures 4.14 and 4.15 are scatter plots of the pole locations for the 
simulated weather plus clutter and clutter only returns for 10 simulated runs. Note 
the existence of weather poles in Figure 4.14 which are not present in the clutter 
only simulations. The decision boundary has been empirically determined based on 
analysis of the pole locations and the apparent clustering of weather poles and non- 
weather poles. As the quantity and the types of data returns increases, a thorough 
analytical derivation of decision boundaries may be undertaken. The boundary in 
Figures 4.14 and 4.15 is a second order polynomial of the form 


, , v 2 8 

g(v ) = 1 

' 1.2u 2 + 450 v 2 + 12 


where v is Doppler velocity and g(v) in the corresponding pole magnitude. The 
windspeed estimate for each range cell can thus be modeled by 


v = 


for 9{vi) < |/*i | 2 

0 else 


(4.2) 


Figure 4.16 details the complete algorithm which derives a windspeed estimate from 
an IQ data series using the two class pattern recognition classification procedure 
presented in this chapter. 

By overlaying the windspeed estimate derived from a second order modal classi- 
fier along with a filtered pulse-pair windspeed estimate onto a bubble plot of the pole 
values, several characteristics of the classifier can be seen. In Figures 4.17 and 4.18 
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Figure 4.1-5 Polar Values of Pole 1 for 10 Simulated Clutter Only Returns 
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Figure 4.16 Two Class Pattern Recognition Windspeed Estimator 
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observe that the windspeed estimate either passes through the pole with the great- 
est Doppler velocity, or it is set to zero. The pulse-pair estimator has a bias toward 
Doppler zero causing its windspeed estimates to be smaller than the modal estimates. 
Also notice that the poles with large Doppler shifts created through pole spreading 
have been removed by the classifier due to their small magnitudes. The pulse-pair 
windspeed estimates also contain some amount of editing from the pulse-pair spectral 
width estimator. When the spectral width becomes larger that a certain threshold 
value, such as 8 m/s in this case, the pulse-pair windspeed estimator is set to zero 
because the spectrum lacks a dominant spectral mode which causes the pulse-pair 
estimator to become erratic. 

4.4 Second Coefficient Windspeed Estimator 
By making the assumption that clutter modes will be primarily near Doppler 
zero, a simplification to the modal windspeed procedure can be achieved. Begin with 
the second order characteristic polynomial from (3.15) and solve for the coefficients 
a-i and a 2 in terms of pole values 


a i — ~(m i + ^2) 

(4.3) 

a 2 = fi! ■ fj, 2 ■ 

(4.4) 

Express (4.4) in polar form 

0-2 = |pi||/u 2 | {cos(0! + $ 2 ) + jsin(/9i + 0 2 )} 

(4.5) 

where 

0i = arg(/ii) 

(4.6) 

6 2 = arg(p; 2 ) . 

(4.7) 


By assuming that one pole will remain near Doppler zero, # 2 , which is defined 
to have the smaller Doppler shift, will be nearly zero. With this assumption the 
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argument of a 2 which is the sum of the arguments of both poles will be nearly equal to 
the argument of pole 1 (pi). The implication of this is that a savings of computational 
resources may be achieved with this approximation as shown in Figure 4.19. For this 
method the decision boundary can be applied directly to the polar value of the second 
AR coefficient a 2 without the need for computing poles or selecting pole 1 and pole 2 
based on Doppler velocities. 
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Figure 4.19 Second Coefficient Windspeed Estimator 
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CHAPTER 5 


RESULTS 

In this chapter windspeed estimates and associated hazard factor values will 
be compared between results from three different windspeed estimators. The pulse- 
pair method will be compared to the two pole classifier and the second coefficient 
windspeed estimator, both described in Chapter 4. The data set used for analysis is 
of a microburst in heavy rainfall which was recorded over Orlando, Florida in 1991. 
This encounter with the microburst has been labeled event 143. From the flight data 
42 frames, or radar returns containing 96 complex IQ samples for each range cell are 
included. The frames were recorded 2.8 seconds apart and all frames were recorded 
when the antenna sweep pattern was aimed directly along the heading of the aircraft. 
Each frame contains 91 range cells numbered 6-96 where range cell 6 is approximately 
0.8 km in front of the aircraft and range cell 96 is approximately 14 km in front of the 
aircraft. The flight began in frame 1 with the aircraft about 8.5 km away from the 
microburst event heading toward the event. At approximately frame 35 the aircraft 
penetrated the center of the microburst outflow. 

Also included in this chapter is an introduction to the fixed range time scan 
method of computing the hazard factor. In this method the radar returns data from 
one range cell which is a constant distance from the aircraft rather than scanning 
throughout a series of range cells. Information is gathered as the aircraft moves 
through a region in space as opposed to a snapshot of multiple range cells which 
return information about a region of space at nearly one moment. Advantages and 
disadvantages of such a method will be discussed as well as a comparison of results 
between the fixed range time scan and the more convention radar snapshot. 



5.1 Comparison of Windspeed Estimates 

Begin the analysis by comparing windspeed estimates between techniques. In 
order to include all 42 frames of data, a three dimensional view of windspeed estimates 
is used. The pulse-pair estimates of windspeed provided the pulse-pair spectral width 
is less than 8 m/s are provided in Figures 5.1 and 5.2. The estimates have been set 
to zero for width values greater than the threshold. Figure 5.1 is a view of actual 
values for each range cell for each frame while in Figure 5.2 the estimates have been 
smoothed by averaging over a three cell by three frame region with the center cell 
and frame reflecting the average. In the smoothed plot windspeed trends are more 
readily seen which will aid in comparing estimation techniques. By observing these 
windspeed trends, characteristics of the microburst event can be seen. First, notice 
that the windspeed features are slanted in a diagonal direction from low frame and 
high range cell to a high frame number and a low range cell. This is created by the 
aircraft motion as it approaches and enters the microburst. At the first snapshot 
(frame 1) the event is seen in the larger range cells because the flight begins at a 
distance of about 8.5 km. During each successive frame the aircraft flies closer to the 
event which shifts its corresponding radar return in to closer range cells. 

Another feature that can be seen in the windspeed plot is a valley that is located 
in front of a rising ridge. The valley is created by negative Doppler shifts which 
indicate a performance increasing headwind and the ridge indicates positive Doppler 
shifts which are induced by a performance decreasing tailwind. By putting these two 
features together while looking at the estimates from one frame only, the characteristic 
“S' 1 -curve pattern appears. These windspeed “slices” for each frame are included in 
Appendix B. 

For a two pole classifier which will be applied to the Orlando event 143. begin 
by observing the distribution of pole 1 for all frames, given in Figure 5.3. Notice that 
the boundary g[v) which was determined for the simulated situation in Chapter 4 
has been overlaid onto the pole distribution for this situation as well. The feature 
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Figure 5.1 Pulse-Pair Velocity Estimates for Event 143 
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Figure 5.2 Pulse-Pair Velocity Estimates for Event 143. smoothed 
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space boundary in this analysis will reflect simulated data on account of the lack of 
any clutter only data available at the present time for the Orlando area. A decision 
boundary determined from a situation which does not have both weather plus clutter 
and clutter only data will most likely not yield results that are as good as when both 
cases are available. For this reason, the classifiers presented in this chapter will use 
the decision boundary determined in Chapter 4. 

As the volume of useful clutter data increases, analyses of pole value distributions 
at certain locations, antenna elevations, and azimuths, may indicate that one decision 
boundary for a particular location is not sufficient. It may be beneficial to implement 
several decision boundaries depending upon these or other factors. At the present 
time it is not known how pole distributions change, given these factors. This analysis 
will be left for future research as more data becomes accessible. 

By implementing the two pole classifier as specified in Chapter 4 with the decision 
boundary as shown in Figure 5.3, the estimated values of windspeed are computed 
and shown in Figures 5.4 and 5.5. Again, Figure 5.5 contains a smoothed sample of 
the windspeed estimates. It is now possible to compare the two pole classifier with the 
pulse-pair estimator over several range cells. Notice that the general shape of both 
estimators is similar in that both contain the diagonal valley and ridge. To compare 
windspeed magnitudes refer to the individual frame plots included in Appendix B. 

One feature that is present in the two pole classifier that is not present in the 
pulse-pair estimator is peaks of large Doppler velocities which occur primarily in the 
range cells following the microburst event. These peaks result from range cells where 
pole spreading occurs and the classifier assigns the spread pole as weather instead of 
recognizing it, as clutter. This error can be seen in the pole scatter plot (Figure 5.3) as 
those poles with a large Doppler velocity and are large in magnitude. These poles do 
not appear to fit with the cluster of weather poles, seen as large magnitude poles with 
Doppler velocities between ±10 m/s. The classifier can not differentiate between these 
poles and weather poles due to them being located in a region where weather poles 
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Figure 5.4 Two Pole Classification Estimates for Event 143 
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Figure 5.5 Two Pole Classification Estimates for Event 143. smoothed 
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may appear. Rejecting these poles by altering the decision boundary will enhance 
estimation for this particular case, but may remove important weather poles when 
used in another situation. 

The next windspeed estimator that will be examined is the second coefficient 
estimator. As shown in Chapter 4 this is a simplified version of a two pole classifier 
which assumes at least one pole will be near Doppler zero for a reliable windspeed 
estimate. In this method the decision boundary is applied to the second AR coeffi- 
cient rather than the pole with the greatest Doppler velocity. Since the second AR 
coefficient is the product of the two complex poles, the difference between the second 
coefficient and pole 1 is the value of pole 2. For returns with mostly stationary clutter 
the value of pole 2 will be approximately 1. Figure 5.6 is a scatter plot of second 
coefficient values which closely matches the scatter plot formed from pole 1 values 
(Figure 5.3). By applying the same decision boundary to the second coefficient that 
was applied to pole 1 in the two pole classifier, the windspeed estimates given in 
Figures 5.7 and 5.8 are produced. 

A comparison between the unsmoothed windspeed estimates from the second co- 
efficient and the two pole classifier reveals that the two estimators yield similar results 
for this situation. Some differences that can be noted are the apparent reduction of 
misclassifications in the second coefficient estimator which results in fewer isolated 
peaks in the velocity estimate plots. This may be attributed to the second coefficient 
having a smaller magnitude than pole 1 as a result of its multiplication by pole 2 
which has a magnitude less than unity. Having a slightly smaller magnitude will 
tend to increase the probability that the second coefficient will be below the decision 
boundary threshold and will be set to zero when pole 1 may be above the threshold. 
It can be observed in the single frame plots included in Appendix B that the second 
coefficient estimator is more often set to zero than the two pole classifier. This rela- 
tion may indicate that separate decision boundaries may be required for each type of 
windspeed estimator. 
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Figure 5.6 Values of Second AR Coefficient for Event 143 







Figure 5.7 Second Coefficient Velocity Estimates for Event 143 



Figure 5.S Second Coefficient Velocity Estimates for Event 143, smoothed 


55 


By further comparing the windspeed estimates provided in Appendix B, several 
more observations can be made about the filtered pulse-pair, the two pole classi- 
fier, and the second coefficient estimator. First, notice that the two modal analysis 
methods generate similar results which indicates that a windspeed estimate computed 
directly from AR coefficients may be feasible. Compared to the pulse-pair method, 
modal analysis tends to generate larger windspeed estimates, which can be attributed 
to the clutter bias present in the filtered pulse-pair estimates. On the other hand, 
the modal analysis results appear to be more erratic than the pulse-pair estimates. 
This is seen more in both the first group of range cells and the last group of range 
cells where the microburst return is not as stong. The erratic behavior of the modal 
estiates may be improved by implementing a more accurate decision boundary and 
forming windspeed estimates from pole values beyond the current range cell. 

5.2 Comparison of Hazard Factors 

The hazard factor values computed in this discussion are derived from the weight- 
ed least-squares procedure discussed in Chapter 2. The hazard factor is computed 
for each frame and the results are recombined to produce another three dimensional 
plot similar to the windspeed estimate plots presented earlier. Figures 5.9 - 5.14 
are included and contain hazard factor values for the windspeed estimators being 
compared. For clarity, only hazard factor values which exceed a threshold of 0.1 are 
included. The remaining values are set to zero in the plot. In this wav, the regions 
with a possible hazard are emphasized. Also presented with the three dimensional 
surface plot of hazard factor is a contour plot of the same data. The contour plot 
assists in locating exactly which frames and corresponding range cells contain a hazard 
factor above the threshold. 

In comparing between the hazard factor values for all three windspeed estimators, 
all reflect a diagonal line where the hazard factor exceeds the threshold. This line 
corresponds to the diagonal valley and ridge seen in all windspeed estimates. The 
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hazard line is a result of a detected headwind (the valley) which is folio ed by a 
tailwind (the ridge). The hazard factor is an attempt to measure the amount of 
performance loss which an aircraft could expect to encounter by proceeding through 
the forward-looking region. 

Upon closer examination of the hazard factor values, some differences between 
the three techniques become apparent. One of the more obvious differences may be 
the existence of two parallel hazard lines in the modal analysis techniques while the 
pulse-pair shows only one diagonal hazard factor line. This can be attributed to tlm 
modal classifier assigning zero windspeed estimates to weather poles in a region near 
Doppler zero. This creates a flat segment in the windspeed estimates where the “S”- 
curve crosses zero velocity. The hazard factor is reduced when it encounters this flat 
segment, but is large on either side where the windspeed estimates diverge. Also, 
the pulse-pair estimator includes fewer spurious peaks in the hazard factor, while 
the two pole classifier includes the most spurious hazard peaks. When comparing 
magnitudes of the hazard factors, all contain values between 0.1 and 0.2 for the 
diagonal lines. A difference does exist, however, in the magnitudes of the isolated 
hazard peaks. The pulse-pair estimator not only has fewer isolated peaks, but these 
peaks have magnitudes that are smaller than the peaks originating from the modal 
analysis methods. 


5.3 Fixed Range Time Scan 

In a typical pulse Doppler radar system the signal processing involved will work 
with the data in the form of energy returns from a number of incremental range 
cells, each a specified distance from the radar antenna. In this way, the processing 
of the data is nearly simultaneous over a spatial region. For radar systems that 
have a moving platform, such as an airborne system, a different technique may be 
employed that analyzes data by fixing the radar scan range to one spatial distance 
and processing additional data over time as the aircraft moves through the spatial 
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region. This type of processing will be referred to as a fixed range time scan. A 
depiction of a both a, snapshot return and a fixed range time scan return are given in 
Figures 5.15 and 5.16. 

One method for understanding how a fixed range time scan functions in com- 
parison to a snapshot is to look again at windspeed estimates for several range cells 
and several frames such as Figures 5.9, 5.11 and 5.13. Recall that a typical snapshot 
would process all range cells one frame at a time. This can be visualized by separat- 
ing the three dimensional plot into a series of two dimensional plots that contain one 
frame number each. The slicing of the plot in this way is illustrated in Figure 5.17. 

The fixed range time scan processes data from a fixed distance over time. This 
can be visualized by slicing the three dimensional plot into individual range cells with 
all frames included. Figure 5.18 is an illustration of this method. The frames become 
a time axis due the sequential order in which they are recorded. For this situation, 
which is a two pole classification estimator of the Orlando event, each range cell 
represents approximately 2.8 seconds, which is the amount of time required for the 
antenna sweep pattern to cross zero azimuth. Because the valley and ridge features of 
the microburst are diagonal to both frame and range cell, it is possible to take slices 
from either a fixed range (cell) or a fixed time (frame). 

One advantage of a fixed range time scan is that all returns will occur from the 
same range which balances the quality of the return across the spatial region. In 
a snapshot, returns from the farther range cells will suffer from signal attenuation 
while the closer range cells contain strong returns. The fixed range time scan allows 
the return data from a spatial region to be processed as soon as it enters a specified 
region of the aircraft. In this way, processing can be compared to an in situ windspeed 
estimator that is positioned at a fixed distance in front of the aircraft. The processing 
of data to derive a hazard factor would be similar to an in situ hazard processor. More 
than one fixed range could be specified in order to allow for different warning intervals 
based on the ability of the radar system to detect a hazard. Strong hazards could be 
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Figure 5.15 Snapshot 



Figure 5.16 Fixed Range Time Scan 
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detected with a greater warning time while events with weaker returns might require 
a closer range and thus, less warning time. 

Finally, the fixed range time scan could work in conjunction with typical radar 
snapshots by providing an additional viewpoint to the current situation. By providing 
feedback to the information gained by radar snapshots, a fixed range time scan could 
increase the reliability of the estimates generated. 
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CHAPTER 6 


CONCLUSIONS 

With conventional signal processing techniques the ability of an airborne pulse 
Doppler radar system to successfully detect the presence of a hazardous windshear 
condition such as a microburst relies heavily upon the correct estimation of windfield 
characteristics. The presence of ground clutter for aircraft operating at low altitudes 
in and near urban terminal areas inhibits the radar’s ability to clearly resolve wind- 
speed estimates. Typically the solution to this problem is to attempt to remove the 
effect of clutter with clutter rejection filters. With Doppler spectrum modal analysis 
as presented and discussed in this thesis, the clutter is not removed by filtering, but 
modeled along with any weather return and discarded through a mode classifier. If 
classified correctly, the weather mode has the potential to yield a windspeed estimate 
which lacks bias from the clutter mode due to the separate modeling of the clutter 
return. 

Modal analysis of Doppler radar data applied to windspeed estimation possesses 
several favorable attributes over more conventional methods such as the pulse-pair 
estimator. Among these attributes are a low computational load due to the absence 
of clutter rejection filtering and the ability to yield unbiased windspeed estimates 
even in low signal-to-clutter environments. Unfavorable attributes of modal analysis 
include variance in the windspeed estimates resulting from misclassification of some 
poles, which increases the probability of false alarm. An improved pole classification 
scheme should be addressed in future research. 

Currently, one disadvantage of modal analysis is the lack of an optimal method 
for determining the decision boundary in the pole classifier. This is primarily due 
to the limited amount of flight data available from which to form a statistical model 
of the poles. As flight data becomes more accessible, more complete knowledge of 



clutter characteristics can be obtained, which should increase the effectiveness of pole 
classifiers. Classification beyond the scope of a single range cell should add stability 
to windspeed estimates as well as reduce the adverse effects of pole spreading and 
discrete clutter. 

The classification nature of modal analysis lends itself to advanced areas such 
as pattern recognition and neural computing. Pattern recognition routines may be 
able to detect hazardous windshear conditions through direct processing of radar IQ 
data without the need for computing windspeed estimates. For example, a pattern 
recognition routine formulated as a data transformation algorithm might be trained 
to recognize the “S”-curve characteristic of a microburst. Then windshear signatures 
across time and/or azimuth can be combined and averaged to improve the ability to 
recognize and detect a hazardous condition. 

It may also be possible to formulate a more robust windspeed estimator based on 
a combination of two or more individual estimators. Because each type of estimator 
has certain strengths and weaknesses, an integration of several estimators could lead 
to better results by using the best aspects of each system. For example, by using 
pulse-pair estimates in conjunction with a two pole model, the confidence in detecting 
hazards may be increased. 

Finally, by examining data presented through a fixed range time scan, windfield 
estimates can be obtained in a method analogous to a predictive in situ hazard de- 
tector. Since each estimate is obtained as the aircraft is moving through a spatial 
region with a constant radar range, signal attenuation remains nearly constant. In 
comparison, a radar snapshot will contain close range cells with little attenuation and 
farther range cells with much attenuation. 
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APPENDICES 



Appendix A 

Autoregressive Modeling 


Autoregressive (AR) modeling is a popular form of spectral estimation due to 
a relatively low computational power required and an improved spectral resolution 
over conventional FFT approaches [16]. AR models are an all-pole model having the 
following spectral density: 


VarU) = 


cr 2 At 

p 

1 + akexp(-j2TrfkAt ) 

k=\ 


(A.l) 


where At is the sampling interval of the process. An AR model of order p requires the 
estimation of the parameters {aj, < 22 , . . . , a p , cr 2 }, where p can be any positive integer. 
Presented here will be the Levinson-Durbin algorithm for fitting the AR parameters 
to a set of sample data. Other methods of AR modeling exist and are discussed by 
Kay and Marple [16] and Keel [9]. 


A.l Levinson-Durbin Algorithm 

This is an iterative procedure for computing the pth order AR parameters asso- 
ciated with a complex series x n based on the relationship between the AR parameters 
and the autocorrelation function of the complex series. In matrix form the relation- 
ship is described by 
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where the autocorrelation estimate of x n is given by 
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This relationship can assume several forms, but in general is known as the Yule- 
Walker equations [16]. The Levinson-Durbin algorithm iterates through successively 
higher order AR parameter sets {an, }, {a 2 i, 0 , 22 , <^ 2)1 ■■■■> { a pi> a p 2 > •• m^} where 
c.ki refers to the ith coefficient of the kth order AR model. The procedure stops when 
the desired model order is reached. Begin the algorithm with the first order model 
parameters 


a 11 


R xx ( 1 ) 

R xx { 0) 

(l-\a n \ 2 )R xx (0) . 


(A. 4) 
(A.5) 


Higher order parameters are recursively computed as follows: 
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The recursive nature of the Levinson-Durbin algorithm has the advantage of com- 
puting not only the AR coefficients for the given model order, but also for all lower 
order models as well, which can help in the determination of the correct model order. 
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Appendix C 

NASA Model Parameters 


SIMULATION PARAMETERS 

Microburst rotation angle 
A/C Distance to Touchdown (km) 
Aircraft Velocity (kts) 
Glideslope Angle (deg) 

No. of Complete Scans 
Time Between Scans (sec.) 

Roll Attitude (deg) 

Pitch Attitude (deg) 

Yaw Attitude (deg) 

Az Integration Range/2 (deg) 

Az Integration Increment (deg) 

El Integration Range/2 (deg) 

El Integration Increment (deg) 
Rng Integration Increment (m) 
Random Number Seed (0-1) 

Runway Number 
Right (1) or Left (2) 

Hazard calc, residual threshold 

MICROBURST & CLUTTER 

Along Track Offset from TD (km) 
Cross Track Offset from TD (km) 
Rain Standard Deviation (m/s) 
Clutter Standard Deviation (m/s) 
Clutter Calc. Flag (l=on,0=off) 
Discrete Calc. Flag (l=on,0=off) 
Reflectivity Calc. Thres. (dBz) 
Minimum Reflectivity (dBz) 

Attenuation Code (0,1,2) 


-90.0 

7.4 

180.0 

3.0 

1.0 
5.0 
0.0 
0.0 
0.0 
5.0 
0.2 

5.0 
0.2 

100.0 

0.224 

26.0 

2.0 
0.1 


- 1.1 

- 0.1 

1.0 

0.5 

1.0 

1.0 

-24.0 

-40.0 (weather & clutter) 
100.0 (clutter only) 

2.0 
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RADAR PARAMETERS 


Initial Radar Range (km) 1.2 

Number of Range Cells 91.0 

Antenna Az - if no scan (deg) -10.0 

Azimuth Scan Range/2 (deg) 0.0 

Azimuth Scan Increment (deg) 10.0 

Antenna Elevation (deg) 1.0 

Transmitted Power (watts) 200.0 

Frequency (GHz) 9.3 

Pulse Width (microsecs) 0.96 

Pulse Interval (microsecs) 266.31 

Receiver Noise Figure (dB) 4.0 

Receiver Losses (db) 3.0 

Antenna Type (1, 2, or 3) 3.0 

Antenna Radius (m) 0.3048 

Aperture Taper Parameter 0.316 

RMS Trans. Phase Jitter (deg) 0.2 

RMS Trans. Freq Jitter (Hz) 0.0 

SIGNAL PROCESSING 

Number of Pulses 96.0 

Number of A/D bits 12.0 

AGC Gain Factor 0.6 

Processing Threshold (dB) 4.0 

Clutter Filter Code (-2 to N) 0.0 

Clutter Filter Cutoff (m/s) 3.0 

No. of Bins for F-f actor Avr. 5.0 
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